System and computer-implemented method for evaluating integrals using stratification by rank-1 lattices

ABSTRACT

A system for numerically evaluating an integral of a function over an s-dimensional integration domain is described, the system comprising a sample point generator, a function value generator and an integral value estimate generator. The sample point generator is configured to generate a selected number of sample points over the integration domain, the sample points being generated such that there is at least one sample point in each of a plurality of strata distributed over the integration domain, the strata being defined by a rank-I lattice. The function value generator is configured to, for respective ones of the sample points, generate a value for the function at the respective sample point. The integral value estimate generator is configured to use the function values generated by the function value generator at the respective sample points in generating an estimate for the value of the integral. The system finds utility in a number of areas, including computer graphics.

FIELD OF THE INVENTION

The invention relates generally to the field of systems andcomputer-implemented methods for evaluating integrals, and moreparticularly to such systems and computer-implemented methods thatevaluate integrals using Monte Carlo and quasi-Monte Carlomethodologies. The invention particularly provides a new and improvedsystem and computer-implemented method for evaluating integrals using aMonte Carlo methodology in which the integration domain is stratifiedusing a stratification methodology in which the number of strata isindependent of the number of dimensions in the integration domain.Specifically, the invention provides a methodology that makes use ofstratification based on rank-1 lattices. Systems andcomputer-implemented methods according to the invention find utility ina number of applications, including but not limited to computergraphics.

BACKGROUND OF THE INVENTION

In computer graphics, a computer is used to generate digital data thatrepresents the projection of surfaces of objects in, for example, athree-dimensional scene, illuminated by one or more light sources, ontoa two-dimensional image plane, to simulate the recording of the imageby, for example, a camera. The simulated camera may include a lens forprojecting the image of the scene onto the image plane, or it maycomprise a pinhole camera in which case no lens is used. Thetwo-dimensional image is in the form of an array of picture elements(which are variably referred to as “pixels” or “pels”), and the digitaldata generated for each pixel represents the color and luminance of thescene as projected onto the image plane at the point of the respectivepixel in the image plane. The surfaces of the objects may have any of anumber of characteristics, including shape, color, specularity, texture,and so forth, which are preferably rendered in the image as closely aspossible, to provide a realistic-looking image.

Generally, the contributions of the light reflected from the variouspoints in the scene to the pixel value representing the color andintensity of a particular pixel are expressed in the form of the one ormore integrals of relatively complicated functions. Since the integralsused in computer graphics generally do not have a closed-form solution,numerical methods must be used to evaluate them to generate the pixelvalue. Typically, a methodology referred to as the “Monte Carlo” methodhas been used in computer graphics to numerically evaluate theintegrals. Generally, in the Monte Carlo method, to evaluate an integral<ƒ>=∫_([0,1[) _(s) ƒ(x)dx   (1)where f(x) is a real function defined on an integration domaincomprising the “s”-dimensional unit cube [0,1)^(s) (that is, ans-dimensional cube each dimension of which includes “zero,” and excludes“one”),a random point set P_(r)=(ξ₁. . . ξ_(r))⊂[0,1)^(s) of “r”independently realized points ξ_(i), i=1, . . . , r, is generated overthe integration domain. The random points ξ_(i) are used as samplepoints for which sample values f(ξ_(i)) are generated for the functionf(x). An estimate {overscore (ƒ)} for the value of the integral(equation (1)) is generated as $\begin{matrix}{{\left\langle f \right\rangle \approx \overset{\_}{f}} = {\frac{1}{r}{\sum\limits_{i = 1}^{r}{f\left( \xi_{i} \right)}}}} & (2)\end{matrix}$As the number “r” of sample points that are used in generating thesample values f(ξ_(i)) increases, the value of the estimate {overscore(ƒ)} will converge toward the actual value of the integral <ƒ>.Generally, the distribution of values of the estimate {overscore (ƒ)}that will be generated for various values of “r,” that is, for variousnumbers of sample points ξ_(i) are distributed around the actual value<ƒ> with a standard deviation σ which can be estimated by$\begin{matrix}{\sigma = \sqrt{\frac{1}{N - 1}\left( {\overset{\_}{f^{2}} - {\overset{\_}{f}}^{2}} \right)}} & (3)\end{matrix}$if the sample points ξ_(i) used to generate the sample values f(ξ_(i))are statistically independent, that is, if the points ξ_(i) are trulypositioned at random in the integration domain. This yields theprobabilistic error bound $\begin{matrix}{{{Prob}\left( {{{{\int_{{\lbrack{0,1})}^{s}}{{f(x)}{\mathbb{d}x}}} - {\frac{1}{r}{\sum\limits_{k = 1}^{r}{f\left( \xi_{k} \right)}}}}} \leq \frac{3{\sigma(f)}}{\sqrt{r}}} \right)} \approx {0.997.}} & (4)\end{matrix}$The rate of convergence of $O\left( \frac{1}{\sqrt{r}} \right)$is independent of the smoothness of the integrand f(x) and the number“s” of dimensions comprising the integration domain.

Several strategies can be used to improve convergence of the estimate{overscore (ƒ)} to the actual value <ƒ> of the integral. One strategy,which is described in the aforementioned Grabenstein patent application,is to use a low-discrepancy, strictly deterministic point methodology togenerate the set of sample points that will be used in generating thesample values f(ξ_(i)) for the estimate. A low-discrepancy, strictlydeterministic methodology will ensure that the sample points aredistributed throughout the integration domain [0,)^(s) without clumpingin particular regions, which can occur if the sample points are randomlydistributed.

Another strategy to improve the convergence is to stratify theintegration domain [0,1)^(s) into a number “K” of disjoint strata A_(j),such that [0,1)²=∪_(j+1) ^(K)A_(j), and separately evaluate theresulting integrals $\begin{matrix}{{\int_{{\lbrack{0,1})}^{s}}{{f(x)}{\mathbb{d}x}}} = {\sum\limits_{j = 1}^{K}{\int_{A_{j}}{{f(x)}{\mathbb{d}x}}}}} & (5)\end{matrix}$using a sample point ξ_(i) obtained from each stratum A_(j). Stratifyingthe integration domain [0,1)^(s) can also be helpful in ensuring thatsample points are reasonably well distributed throughout the integrationdomain [0,1)^(s), and that they do not clump in the integration domain,which can occur if the locations of the sample points ξ_(i) are randomlylocated in the integration domain.

Several stratification strategies have been proposed. Typically,previously-proposed stratification strategies, particularlystratification strategies that are based on axis-aligned grids, sufferexponential growth in the number of strata with increasing numbers ofdimensions comprising the integration domain. Since integrands ƒ thatare used in fields such as, for example, computer graphics, often havevery high-dimensional integration domains, stratifying domains usingsuch stratification strategies can be computationally quite intensive.

SUMMARY OF THE INVENTION

The invention provides a new and improved system andcomputer-implemented method for evaluating integrals using a Monte Carlomethodology in which the integration domain is stratified using astratification methodology in which the number of strata is independentof the number of dimensions in the integration domain. Specifically, theinvention provides a methodology that makes use of stratification basedon rank-1 lattices.

In brief summary, in one aspect, the invention provides a system fornumerically evaluating an integral of a function over an s-dimensionalintegration domain, the system comprising a sample point generator, afunction value generator and an integral value estimate generator. Thesample point generator is configured to generate a selected number ofsample points over the integration domain, the sample points beinggenerated such that there is at least one sample point in each of aplurality of strata distributed over the integration domain, the stratabeing defined by a rank-1 lattice. The function value generator isconfigured to, for respective ones of the sample points, generate avalue for the function at the respective sample point. The integralvalue estimate generator is configured to use the function valuesgenerated by the function value generator at the respective samplepoints in generating an estimate for the value of the integral.

In another aspect, the invention provides a computer implemented methodof numerically evaluating an integral of a function over ans-dimensional integration domain, the method comprising a sample pointgeneration step, a function value generation step and an integral valueestimate generation step. During the sample point generation step, aselected number of sample points are generated over the integrationdomain, the sample points being generated such that there is at leastone sample point in each of a plurality of strata distributed over theintegration domain, the strata being defined by a rank-1 lattice. Duringthe function value generation step, for respective ones of the samplepoints, a value for the function at the respective sample point isgenerated. During the integral value estimate generation step, thefunction values generated during the function value generation step atthe respective sample points are used in generating an estimate for thevalue of the integral.

In another aspect, the invention provides computer program product foruse with a computer to provide a system for numerically evaluating anintegral of a function over an s-dimensional integration domain. Thecomputer program product comprises a computer readable medium havingencoded thereon a sample point generator module, a function valuegenerator module and an integral value estimate generator module. Thesample point generator module is configured to enable the computer togenerate a selected number of sample points over the integration domain,the sample points being generated such that there is at least one samplepoint in each of a plurality of strata distributed over the integrationdomain, the strata being defined by a rank-1 lattice. The function valuegenerator module is configured to enable the computer to, for respectiveones of the sample points, generate a value for the function at therespective sample point. The integral value estimate generator module isconfigured to enable the computer to use the function values generatedby the function value generator at the respective sample points ingenerating an estimate for the value of the integral.

BRIEF DESCRIPTION OF THE DRAWINGS

This invention is pointed out with particularity in the appended claims.The above and further advantages of this invention may be betterunderstood by referring to the following description taken inconjunction with the accompanying drawings, in which:

FIG. 1 depicts an illustrative computer graphics system that evaluatesintegrals using a Monte Carlo methodology in which the integrationdomain is stratified using a stratification methodology in which thenumber of strata is independent of the number of dimensions in theintegration domain;

FIG. 2 depicts a table of information that is useful in understandingthe operation of the computer graphics system depicted in FIG. 1.

DETAILED DESCRIPTION OF AN ILLUSTRATIVE EMBODIMENT

The invention provides a computer graphics system and method forgenerating pixel values for pixels in an image of a simulated scene thatmakes use of a Monte Carlo methodology in which an integral or integralsis evaluated of function(s) by using sample points determined bystratifying the integration domain using rank-1 lattices. Thefunction(s) represent the contributions of the light that is emittedfrom simulated light sources and reflected from the various simulatedsurfaces in the scene directed towards a simulated camera, and theintegral(s) provide the pixel values for the respective pixels in theimage. Stratifying the integration domain in this manner provides thatthe number of strata will not grow exponentially with the number ofdimensions in the integration domain, which can occur with, for example,axis-aligned strata.

FIG. 1 attached hereto depicts an illustrative computer system 10 thatmakes use of such a stratification methodology. With reference to FIG.1, the computer system 10 in one embodiment includes a processor module11 and operator interface elements comprising operator input componentssuch as a keyboard 12A and/or a mouse 12B (generally identified asoperator input element(s) 12) and an operator output element such as avideo display device 13. The illustrative computer system 10 is of theconventional stored-program computer architecture. The processor module11 includes, for example, one or more processor, memory and mass storagedevices, such as disk and/or tape storage elements (not separatelyshown), which perform processing and storage operations in connectionwith digital data provided thereto. The operator input element(s) 12 areprovided to permit an operator to input information for processing. Thevideo display device 13 is provided to display output informationgenerated by the processor module 11 on a screen 14 to the operator,including data that the operator may input for processing, informationthat the operator may input to control processing, as well asinformation generated during processing. The processor module 11generates information for display by the video display device 13 using aso-called “graphical user interface” (“GUI”), in which information forvarious applications programs is displayed using various “windows.”Although the computer system 10 is shown as comprising particularcomponents, such as the keyboard 12A and mouse 12B for receiving inputinformation from an operator, and a video display device 13 fordisplaying output information to the operator, it will be appreciatedthat the computer system 10 may include a variety of components inaddition to or instead of those depicted in FIG. 1.

In addition, the processor module 11 includes one or more network ports,generally identified by reference numeral 14, which are connected tocommunication links which connect the computer system 10 in a computernetwork. The network ports enable the computer system 10 to transmitinformation to, and receive information from, other computer systems andother devices in the network. In a typical network organized accordingto, for example, the client-server paradigm, certain computer systems inthe network are designated as servers, which store data and programs(generally, “information”) for processing by the other, client computersystems, thereby to enable the client computer systems to convenientlyshare the information. A client computer system which needs access toinformation maintained by a particular server will enable the server todownload the information to it over the network. After processing thedata, the client computer system may also return the processed data tothe server for storage. In addition to computer systems (including theabove-described servers and clients), a network may also include, forexample, printers and facsimile devices, digital audio or video storageand distribution devices, and the like, which may be shared among thevarious computer systems connected in the network. The communicationlinks interconnecting the computer systems in the network may, as isconventional, comprise any convenient information-carrying medium,including wires, optical fibers or other media for carrying signalsamong the computer systems. Computer systems transfer information overthe network by means of messages transferred over the communicationlinks, with each message including information and an identifieridentifying the device to receive the message.

It will be helpful to initially provide some background on operationsperformed by the computer graphics system 10 in generating an image.Generally, the computer graphics system 10 generates an image thatattempts to simulate an image of a three-dimensional scene as would berecorded by a camera. In that operation, for each pixel of the camera'simage plane, the computer graphics system 10 numerically evaluates oneor more integrals of functions that represent light that is directedtoward the respective pixel from the scene. An illustrative functionwill be described below. The functions typically contain a number ofterms, including one or more terms that represent light that is directlyincident on the respective pixel from one or more light source sources,as well as one or more terms that represent light that, after havingbeen emitted by a light source, is reflected by one or more surfaces ofobjects in the scene and ultimately directed toward the pixel.

To numerically evaluate the integral, the computer graphics system 10shoots, or traces, simulated rays from a plurality of points in thepixel into the scene. At points in the scene at which the shot raysintersect surfaces of objects in the scene, the computer graphics systemdetermines the contributions of light that are directed from thoseintersection points in the scene toward the respective points in thepixel from which the rays were shot. The computer graphics system candetermine the contributions that are directed to the respective pointsin the pixel in a number of ways. In a path tracing methodology, whichis described in more detail below, the computer graphics system 10determines the illumination at the point in the scene at which the rayshot from the pixel intersects an object in the scene by shooting one ormore rays from the intersection point in a plurality of directions, anddetermining the extent to which that point is illuminated directly bylight sources, and indirectly by light that is reflected off of othersurfaces of objects in the scene. After the computer graphics system hasdetermined the extent to which the intersection point is illuminated,the direction or directions from which the point is illuminated, andother information that will be known by those skilled in the art, itgenerates a value that represents the light reflected toward the pointin the pixel from which the ray was shot. The computer graphics systemwill repeat these operations for the various points in the pixel fromwhich rays are shot, and will determine a pixel value for the entirepixel from the various values that are generated therefor. The computergraphics system also repeats the operations for all of the pixels togenerate the final image. In a bidirectional path tracing methodology,the computer graphics system also shoots rays from light source orsources into the scenes and determines the pixel values in relation topoints in at which those rays intersect points on surfaces of objects inthe scene.

In a photon map methodology, prior to shooting rays from the pixels inthe image plane into the scene, the computer graphics system 10 firstgenerates a photon map representing photons that are incident on thesurfaces of objects in the scene directly from the light sources andindirectly by photons that are reflected from various surfaces in thescene. After the computer graphics system has generated the photon map,it will shoot rays from the points on the pixels into the scene, andgenerate the pixel value in relation to the photons, if any, that areproximate the point of intersection of those rays with points onsurfaces of the objects in the scene.

In both methodologies, as well as other methodologies as will beappreciated by those skilled in the art, the computer graphics system 10uses the coordinates of the sample points that are obtained from theintegration domain of the integral to control various aspects regardingthe geometry while generating the image. These aspects may include, forexample, determining the locations of the points on the various pixelsfrom which rays are shot, determining the locations on light sourcesfrom which rays representing photons are shot, determining directions ofreflection, determining whether to terminate a ray at a surface or toallow the ray to be reflected, and other aspects as will be aware tothose skilled in the art.

As noted above, formally, the operations performed by the computergraphics system in generating pixel values for an image comprise thenumerical evaluation of an integral over an integration domain that isdefined over an s-dimensional unit cube [0,1)^(s). In numericallyevaluating the integral, the computer graphics system partitions theintegration domain [0,1)^(s) into a number “K” of disjoint strata A_(j),such that [0,1)^(s)∪=_(j=1) ^(K)A_(j), and separately evaluates theresulting integrals $\begin{matrix}{{\int_{{\lbrack{0,1})}^{s}}{{f(x)}{\mathbb{d}x}}} = {\sum\limits_{j = 1}^{K}{\int_{A_{j}}{{f(x)}{\mathbb{d}x}}}}} & (6)\end{matrix}$using a sample point obtained from each stratum A_(j). The strata A₄ aredefined by an “s” dimensional lattice. The lattice divides theintegration domain into a plurality of strata A_(j) of equal size, andthe computer graphics system 10 obtains a sample point from each stratumA_(j) for use in generating the approximation for the integral.

A lattice is a discrete set of points L:=P_(n)+Z^(s)⊂R^(s) that isclosed under addition and subtraction, where P_(n):={z₀, . . . ,z_(n-1)l }, “Z” is an integer, “R” is a real number, and superscript “s”corresponds to the dimension of the lattice. In the case of a latticedefined by an axis-aligned regular grid, the dimension “s” is the numberof basis vectors that, in linear combinations, yield the lattice points.According to another methodology, and in accordance with the invention,the lattice is defined using a single generator vector g ε N^(s), where“N” is a member of the natural numbers, with the coordinates of thelattice points being defined by $\begin{matrix}{{z_{j}:={\frac{j}{n}g}},} & (7)\end{matrix}$where j=0, . . . ,n-1. The lattice is referred to as a “rank-1” lattice,because it is generated using a single generator vector g.

An example of a rank-1 lattice that will be used herein is a latticethat is generated based on the Fibonacci sequence. The Fibonaccisequence {F_(k)} is defined as F_(k):=F_(k−1)+F_(k−2), with F₁=F₂=1 andindex “k” an integer. The two-dimensional Fibonacci lattice, which isbased on the Fibonacci sequence, for n=F_(k), k≧3, lattice vertices hasthe generator vectorg=(1, F _(k−1))   (8)It will be apparent from equations (7) and (8) that the coordinates(α,β) of the vertices z_(j) comprising a two-dimensional Fibonaccilattice on the unit square [0,1)² can be readily generated in a“for”-loop, in which the first component “α” is $``{j\frac{1}{F_{k}}}"$and the second component “β” is${``{j\frac{F_{k - 1}}{F_{k}}{mod}\quad 1}"},$ since, as noted above,n=F_(k), the “k-th” element of the Fibonacci sequence.Fibonacci lattices of dimension “s” higher than “two” can be defined ina similar manner using an “s”-dimensional generator vector. For example,the generator vectorg = (1, F_(k − 1), F_(k − 1)², …  , F_(k − 1)^(s − 1))can be used, which would provide a lattice in Korobov form.

As noted above, the computer graphics system 10 obtains sample pointsfrom strata A_(j) that are defined by the lattice. The strata relate tounit cells that are associated with the respective vertices z_(j). Theunit cell of a rank-1 lattice is not uniquely specified, but a suitableunit cell, which is referred to as a Delaunay unit cell, can bespecified such that the volume of an s-dimensional sphere that isinscribed in the unit cell is maximized. The unit cells that result onthe rank-1 lattice are all of similar shape and can be mapped from oneto another by means of rigid body transformations that basically performa translation and/or rotation to move them from vertex to vertex.

Each unit cell is identified by the Voronoi diagram of the rank-1lattice, whose nerve is constructed out of the “s” shortest vectors v₁,. . . , v_(s) from the origin of the coordinate system defining thelattice, to points of the lattice, where the value of the firstcomponent of each vector v_(j) ⁽¹⁾ is greater than zero. This impliesthat, for all vectors, v_(j), v_(j)≠0, and further that the value of“n,” the number of vertices of the lattice, is greater than “one.” Thevectors span an s-dimensional parallelopipedB:=(v₁ . . . v_(s))   (9)whose volume ${B}:={{{\det\quad B}} = {\frac{1}{n}.}}$Each of the “n” s-dimensional parallelopipeds in a rank-1 latticerestricted to the unit cube [0,1)^(s) is a stratum A_(j) that is inducedby the bijection $\begin{matrix}{{{R_{j}:\left\lbrack {0,1} \right)^{s}}->A_{j}}\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {B\quad x}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.} & (10)\end{matrix}$where the j-th lattice vertex $z_{j} = {\frac{j}{n} \cdot g}$is used as an offset into the “s”-dimensional unit cube [0,1)^(s) and“B,” which is used as a matrix with the vectors v_(j) comprising thecolumns of the matrix, defines a linear transformation that represents arotation and/or shear. Equation (10) defines what can be referred to asa “reduced Cranley-Patterson rotation,” which determines the location ofthe “j-th” sample point at which the function ƒ will be evaluated forthe contribution to the estimate {overscore (ƒ)} of the value of theintegral. Generally, the first term in the second line of equation (10),namely, ${``{\frac{j}{n} \cdot g}"},$j=0, . . . ,n-1, identifies the locations of the vertices for the basic“unrotated” rank-1 lattice. The second term in the second line ofequation (9), namely, B x, essentially provides a displacement fromlocations of the respective unrotated vertices to the particular samplepoint that will be used in evaluating the function. As will be describedbelow, the computer graphics system 10 may use the same value for “x”for some or all of the vertices (reference the discussion belowregarding correlated sampling), or it may use different values(reference the discussion below regarding jittered sampling). Inequation (10), if different values are to be used for “x,” “x” may bej-th element of an s-dimensional point set P_(r)=(ξ₁, . . .ξ_(r))⊂[0,1)^(s), which may be generated using a random methodology, astrictly-deterministic low-discrepancy methodology, or othermethodologies that will be apparent to those skilled in the art.

The particular parallelopiped that is defined by B will depend on theparticular form of the generator vector g that is used to define thelattice, the number “n” of vertices that are to comprise the lattice,the number of dimensions, and other criteria as will be appreciated bythose skilled in the art. For the two-dimensional Fibonacci lattice withn=F_(k)>1 vertices, the matrix B is given by $\begin{matrix}{B = {\frac{1}{F_{k}}\begin{pmatrix}j_{1} & j_{2} \\{j_{1}F_{k - 1}} & {{- j_{2}}F_{k - 1}}\end{pmatrix}}} & (11)\end{matrix}$where$j_{1} = {{F_{{2 \cdot {\lfloor\frac{k - 1}{4}\rfloor}} + 1}\quad{and}\quad j_{2}} = {F_{2 \cdot {\lfloor\frac{k + 1}{4}\rfloor}}.}}$FIG. 2 depicts a table that, for two-dimensional Fibonacci lattices ofvarious numbers “n” of vertices, gives the generator vector g and valuesof j₁ and j₂ for use in connection with equation (11). As noted above,in the Fibonacci sequence, F₁=F₂=1. Accordingly, the third element ofthe sequence, F₃, is the sum of the two preceding elements, and soF₃=F₂+F₁=2, which, in turn, corresponds to the number “n” of vertices inthe lattice, as shown in the first line of the table depicted in FIG. 2.For the two-dimensional lattice associated with the F₃ element of theFibonacci sequence

-   -   (i) the generator vector g=(1,F_(k−1))=(1,F₂)=(1,1);        ${{{({ii})\quad j_{1}} = {F_{{2 \cdot {\lfloor\frac{3 - 1}{4}\rfloor}} + 1} = {F_{{2 \cdot 0} + 1} = {F_{1} = 1}}}};}\quad$        ${{({iii})\quad j_{2}} = {F_{2 \cdot {\lfloor\frac{3 + 1}{4}\rfloor}} = {F_{2 \cdot 1} = {F_{2} = 1}}}};$        ${{{and}\quad{{so}({iv})}\quad B} = {\frac{1}{2}{\begin{pmatrix}        1 & 1 \\        1 & {- 1}        \end{pmatrix}.}}}\quad$        The values for items (i) through (iii) are depicted in the        second, third and fourth columns, respectively, of the table        depicted in FIG. 2. It will be apparent that the term        $\frac{j}{n} \cdot g$        in the second line of equation (9), which defines the unrotated        lattice, corresponds to vertices having coordinates (0,0) for        j=0 and (1/2,1/2) for j=1. The matrix B, in turn, defines strata        that are on the order of 1/2-by-1/2 in size. That is,

(a) the stratum that is associated with the vertex having coordinates(0,0) comprises the quadrilateral bounded by line segments betweenpoints (0,0), (1/2,0), (1/2,1/2) and (0,1/2), with the quadrilateralincluding the line segments between points (0,0)-(1/2,0) and(0,0)-(0,1/2) (except for points (1/2,0) and (0,1/2)), and not includingthe line segments between points (1/2,0)-(1/2,1/2) and(0,1/2)-(1/2,1/2), and

(b) the stratum that is associated with the vertex having coordinates(1/2,1/2) comprises the quadrilateral bounded by line segments betweenpoints (1/2,1/2), (1,1/2), (1,1) and (1/2,1), with the quadrilateralincluding the line segments between points (1/2,1/2)-(1,1/2) and(1/2,1/2)-(1/2,1) (except for points (1,1/2) and (1/2,1)), and notincluding the line segments through points (1,0)-(1,1) and (0,1)-(1,1).

As noted above, the strata associated with the respective verticesessentially represent a rigid translation, and it is apparent that thestratum associated with the vertex having coordinates (1/2,1/2) has thesame size and shape as the stratum associated with the vertex havingcoordinates (0,0). The particular location within each stratum at whichthe vertex of the rotated lattice will be determined by the value of “x”(reference the second line of equation (10)).

As another illustration of the generation of values for g, j₁, j₁, j₂and B, for the two-dimensional 19 lattice associated with the ninthelement of the Fibonacci sequence, F₉, there are F₉=F₈+F₇=21+13=34vertices and, for the lattice that is associated with th

-   -   (i) the vector generator g=(1,F_(k−1))=(1,F₈)=(1,21);        ${{{({ii})\quad j_{1}} = {F_{{2 \cdot {\lfloor\frac{9 - 1}{4}\rfloor}} + 1} = {F_{{2 \cdot {\lfloor\frac{8}{4}\rfloor}} + 1} = {F_{{2 \cdot 2} + 1} = {F_{5} = 5}}}}};}\quad$        ${{{({iii})\quad j_{2}} = {F_{2 \cdot {\lfloor\frac{9 + 1}{4}\rfloor}} = {F_{2 \cdot {\lfloor\frac{10}{4}\rfloor}} = {F_{2 \cdot 2} = {F_{4} = 3}}}}};{{{and}\quad{{so}({iv})}\quad B} = {{\frac{1}{34}\begin{pmatrix}        5 & 3 \\        {5 \cdot 21} & {{- 3} \cdot 21}        \end{pmatrix}} = {\frac{1}{34}{\begin{pmatrix}        5 & 3 \\        105 & {- 63}        \end{pmatrix}.}}}}}\quad$        It will be apparent that the term $\frac{j}{n} \cdot g$        in the second line of equation (10), which defines the unrotated        lattice, corresponds to vertices having coordinates (0,0) for        j=0, (1/34,21/34) for j=1,(2/34,8/34) for j=2, (3/34,29/34) for        j=3, (4/34,16/34) for j=4, (5/34,3/34) for j=5, and so forth.        The matrix B, in turn, defines for the vertex associated with        j=0, namely, the vertex having coordinates (0,0), the stratum        comprising the quadrilateral bounded by line segments between        points (0,0), (5/34,3/34), (8/34,8/34) and (3/34,5/34), with the        quadrilateral including the line segments between points        (0,0)-(5/34,3/34) and (0,0)-(3/34,5/34) (except for points        (5/34,3/34) and (3/34,5/34), and not including the line segments        between points (5/34,3/34)-(8/34,8/34) and        (3/34,5/34)-(8134,8/34). The strata that are associated with the        other vertices can be determined in a similar manner. As with        the strata associated with the two-vertex lattice (F₃) described        above, in the case of the thirty-four vertex (F₉) lattice, the        strata associated with all of the vertices have the same sizes        and shapes

Other lines of the table depicted in FIG. 2 depict similar informationfor two-dimensional Fibonacci lattices having different numbers ofvertices.

Benefits of stratified sampling using rank-one lattices will beillustrated using several examples. Generally, in contrast tostratification using grid-aligned lattices, stratifying an“s”-dimensional domain using rank-1 lattices provides that the number ofstrata will not grow exponentially with the number of dimensions in theintegration domain. This will be illustrated in connection with jitteredsampling, correlating sampling and correlated trajectory samplingmethodologies.

Jittered Sampling by Rank-1 Lattices

Generally, jittered sampling is achieved by placing the samples pointsused for evaluating the integrand ƒ in random positions withinrespective strata. Using the stratification that is induced by rank-1lattices (equation (10)), the value of an integral can be estimated asfollows: $\begin{matrix}\begin{matrix}{{\int_{{\lbrack{0,1})}^{s}}{{f(x)}{\mathbb{d}x}}} = {\sum\limits_{j = 0}^{n - 1}{\int_{A_{j}}{{f(x)}{\mathbb{d}x}}}}} \\{= {\frac{1}{n}{\sum\limits_{j = 0}^{n - 1}{\int_{{\lbrack{0,1})}^{s}}{{f\left( {R_{j}(x)} \right)}{\mathbb{d}x}}}}}} \\{\approx {\frac{1}{n}{\sum\limits_{j = 0}^{n - 1}{f\left( {R_{j}\left( \xi_{j} \right)} \right)}}}}\end{matrix} & (12)\end{matrix}$In equation (12), the first line reflects the fact that the integralover the entire integration domain, represented by the left hand side,corresponds to the sum of the integrals over the various strata A_(j),j=0, . . . ,n-1. The second line of equation (12) follows from the factthat, since, if “x” is a continuous variable spanning the integrationdomain, R_(j)(x) will also be continuous and span the same entireintegration domain. In that case, if the integral is summed “n” times,the value of the integral can be obtained by dividing the sum by “n.”Finally, the third line reflects that the value of the “j-th” integralin the sum in the second line can be estimated by evaluating thefunction ƒ using the “j-th” sample point, which corresponds toR_(j)(ξ_(j)).

In determining the value for the estimate {overscore (ƒ)}, whichcorresponds to the last line of equation (12), the computer graphicssystem 10 will, for each value of index “j,” make use of the generatorvector g to generate the $\frac{j}{n} \cdot g$term of R_(j)(ξ_(j)) (reference equation (10)), and the “j-th” elementξ_(j) of an “s”-dimensional random sequence {ξ_(j)}_(j = 0)^(n − 1)to generate the Bx term, where the “j-th” element ξ_(j) is used asvector “x” in generating R_(j)(ξ_(j)) (equation (10)). The sum Rj(ξ_(j))is used as the sample point in evaluating the function ƒ for therespective “j-th” sample point. Since different elements ξ_(j) of therandom sequence {ξ_(j)}_(j = 0)^(n − 1)are used in generating R_(j)(ξ_(j)) for different values of index “j,”the sample points ξ_(j) are likely to have different displacementswithin the various strata A_(j). However since there is at least onesample point R_(j)(ξ_(j)) in each stratum, and since the strata areuniformly distributed over the integration domain, the sample pointswill be reasonably well distributed around the integration domain. Itwill be appreciated that several elements ξ_(j) of the random sequence{ξ_(j)}_(j = 0)^(n − 1)may clump, and may even have the same value, but the stratificationinduced by the $\frac{j}{n} \cdot g$term of equation (10) will ensure that the sample points R_(j)(ξ_(j))that are used in evaluating the function ƒ are well distributed over theintegration domain. A more uniform distribution may be achieved if theelements ξ_(j) of the random sequence {ξ_(j)}_(j = 0)^(n − 1)are generated using, for example, a strictly-deterministic,low-discrepancy methodology.

Using rank-1 lattices, the computer graphics system 10 can readilypartition the s-dimensional unit cube [0,1)^(s) into a number “n” ofstrata, with the number “n” being independent of the number “s” ofdimensions. The estimator for the value of the integral, which isdescribed in equation (12), and particularly in the last line of thatequation, is unbiased and can also be used in connection with generatingestimates of the values of integrals of functions related to fieldsother than computer graphics.

Correlated Sampling by Rank-1 Lattices

As with jittered sampling, correlated sampling also involves evaluatinga sum of integrals of a function ƒ over the s-dimensional unit cube.However, unlike jittered sampling, in which elements of a random“s”-dimensional sequence {ξ_(j)}_(j = 0)^(n − 1)is used in generating the sample points at which the function ƒ isevaluated (equation (12)), in correlated sampling one randomly generatedelement ξ is used in generating all of the sample points (equation (10))that are used in evaluating the function ƒ for use in the sum. This isderived from equation (12) by, with reference to the second line ofequation (12), interchanging the finite sum and the integral. That is,$\begin{matrix}\begin{matrix}{{\frac{1}{n}{\sum\limits_{j = 0}^{n - 1}{\int_{{\lbrack{0,1})}^{s}}{{f\left( {R_{j}(x)} \right)}{\mathbb{d}x}}}}} = {\frac{1}{n}{\int_{{\lbrack{0,1})}^{s}}{\sum\limits_{j = 0}^{n - 1}{{f\left( {R_{j}(x)} \right)}{\mathbb{d}x}}}}}} \\{\approx {\frac{1}{n}{\sum\limits_{j = 0}^{n - 1}{f\left( {R_{j}(\xi)} \right)}}}}\end{matrix} & (13)\end{matrix}$It will be appreciated that, in equation (13), the left hand side of thefirst line of the equation corresponds to the right hand side of thesecond line of equation (12). In generating the estimate {overscore(ƒ)}, which corresponds to the last line of equation (13), the computergraphics system 10 will, for each value of “j,” make use of thegenerator vector g in generating the $\frac{j}{n} \cdot g$term of R_(j)(ξ), and the same randomly-generated vector “ξ” as “x” inthe Bx term of the sample point R_(j)(ξ). As noted above, in generatingthe image, the computer graphics system 10 will separately evaluate theintegral for each of the pixels in the image; and in that case thecomputer graphics system 10 may use the same element ξ for all of theintegrals associated with the various pixels, or, alternatively, thecomputer graphics system 10 may use different random elements . . .ξ_(i−1), ξ_(i), ξ₁₊₁, . . for the integral associated with the . . .“i-1st,” “i-th”, “i+1st” pixel, respectively. One advantage of usingdifferent elements for the various integrals is that that may result inthe generation of noise in the image, which can obscure aliasingartifacts that might otherwise be present in an image.

Since stratified samples are independent, there is no correlationbetween them; in particular, there is no minimum distance property.Therefore, no variance reduction will be guaranteed. Minimum distanceconstraints can be imposed by, for example, applying Lloyd's relaxationalgorithm in a conventional manner. Combined with modified sampleweights by post-stratification, the estimate {overscore (ƒ)}converges tothe actual value <ƒ> of the integral. For Riemann-integrable functionsthe unbiased estimator (reference the second line of equation (13)) hasa reduced variance of order${O\left( {\frac{1}{n^{2}}\ln^{2s}n} \right)},$if the set of sample points (R_(j)(ξ))_(j = 0)^(n − 1)is of sufficiently low discrepancy.

The estimator (reference the second line of equation (13)) that is usedfor an integral can be realized by generating one random s-dimensionalvector ξ ε [0,1)^(s) and replicating it using equation (10), that is,using the product of B and the random element ξ as an offset for each ofthe unrotated lattice vertices, which, as noted above, are defined$\frac{j}{n} \cdot {g.}$The resulting correlated samples have a guaranteed minimum distance andgood uniformity properties, if the generator vector g is suitablychosen. If values for “x” are uniformly distributed the transformation Bcan be omitted, in which case the bijection (reference equation (10))reduces to $\begin{matrix}{{{R_{j}^{CP}:\left\lbrack {0,1} \right)^{s}}->\left\lbrack {0,1} \right)^{s}}\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + x} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.} & (14)\end{matrix}$which corresponds to a Cranley-Patterson rotation on a rank-1 lattice.

Correlated Trajectory Splitting by Rank-1 Lattices

Trajectory splitting can improve efficiency depending on the correlationof the integrand with respect to its lower dimensional projections.Trajectory splitting is used in a number of types of simulations,including but not limited to computer graphics. Trajectory splitting canbe used in several ways in computer graphics, one of which will bedescribed here and another below. One use of trajectory splitting inconnection with computer graphics is tracing multiple shadow rays to anlight source that has a non-zero area. As noted above, in ray tracing,path tracing and similar methodologies, to determine a pixel value for apixel in an image, the computer graphics system 10 traces a plurality ofrays from respective points in the pixel into a simulated scene. Foreach such ray, if the ray intersects a surface of an object in thescene, the computer graphics system 10 can trace one or more rays, whichare referred to as shadow rays, from the point of intersection of theshot ray to various points on the area of the light source. If a shadowray intersects another object in the scene, the point of intersection ofthe shot ray is in the other object's shadow for that shadow ray.

On the other hand, if the shadow ray does not intersect another objectin the scene, the light source is deemed to contribute to the directillumination of the point in the scene for the pixel from which the raywas shot. The shadow rays that are traced from the point of intersectionto the light source are distributed over the surface of the lightsource, and so some of the shadow rays may intersect other objects inthe scene, whereas others of the shadow rays may not. The extent towhich the point of intersection of the shot ray is illuminated by thelight source is a function of the number of shadow rays from that pointon the object that do not intersect other objects. Generally, if a scenecontains multiple light sources, similar operations will be performedfor each light source to determine the extent to which the point atwhich the ray that was shot from the pixel intersects the object, isilluminated by the various light sources.

As noted above, a pixel also has an area, and the computer graphicssystem will typically shoot simulated rays into the scene from severalpoints within the area of the pixel, and trace shadow rays from pointsat which they intersect the various objects. Generally, the variance ofthe integrand with respect to the area of the pixel is less than thevariance with respect to the area of the light source.

In trajectory splitting, the samples (illustratively, the shadow rays)that are split off of one instance (illustratively, the ray that istraced from the pixel into the scene) do not need to be independent.Continuing with the computer graphics shadow ray example, using a vector“x” to refer to a location in the pixel, and a vector “y” to refer to alocation on the light source, the function ƒ in the integral is afunction of both vector “x” and vector “y,” that is f(x,y), with vector“x” being associated with s₁ dimensions and vector “y” being associatedwith s₂ dimensions. Accordingly, the portion of the integration domainthat is associated with vector “x” is the s₁-dimensional unit cube[0,1)^(s) ¹ and the portion of the integration domain that is associatedwith the vector “y” is the s₂ dimensional unit cube [0,1)₂ ¹ . In thatcase, if s₁+s₂=s, which will be the case if the vectors “x” and “y” areassociated with all of the dimensions comprising the integration domain,the integral of ƒ is evaluated as follows $\begin{matrix}\begin{matrix}{{\int_{{\lbrack{0,1})}^{s_{1}}}{\int_{{\lbrack{0,1})}^{s_{2}}}{{f\left( {x,y} \right)}{\mathbb{d}y}{\mathbb{d}x}}}} = {\int_{{\lbrack{0,1})}^{s_{1}}}{\sum\limits_{j = 0}^{n - 1}{\int_{A_{j}}{{f\left( {x,y} \right)}{\mathbb{d}y}\quad{\mathbb{d}x}}}}}} \\{= {\int_{{\lbrack{0,1})}^{s_{1}}}{\frac{1}{n}{\int_{{\lbrack{0,1})}^{s_{2}}}{\sum\limits_{j = 0}^{n - 1}{{f\left( {x,{R_{j}(y)}} \right)}{\mathbb{d}y}\quad{\mathbb{d}x}}}}}}} \\{\approx {\frac{1}{r}{\sum\limits_{i = 0}^{r - 1}{\frac{1}{n}{\sum\limits_{j = 0}^{n - 1}{f\left( {\xi_{i},{R_{j}\left( \zeta_{i} \right)}} \right)}}}}}}\end{matrix} & (15)\end{matrix}$where ξ_(i) and ζ_(i) are the “i-th” elements of two random sequences,namely, an s₁-dimensional random sequence {ξ_(i)}_(i = 0)^(r − 1)and an s₂-dimensional random sequence {ζ_(i)}_(i = 0)^(r − 1).Equation (15) essentially reflects the fact that the rank-1 lattice,instead of being used to stratify the entire “s”-dimensional integrationdomain, can be used to stratify a subset s₂<s of the dimensionscomprising the integration domain. Continuing with the shadow-rayexample, the element ξ_(i) is used to determine the location of a pointwithin the pixel from which the simulated ray is to be shot into thescene. On the other hand, R_(j)(ζ_(i)) (reference equation (10)) is usedto determine the location on the area light source toward which a shadowray will be traced from the point at which the ray shot from the pixelintersects an object in the scene.

It will be apparent that, in evaluating the integral according to thelast line of equation (15), the sampling methodology essentiallycorresponds to subdividing the area light source into a number j=0, . .. ,n-1 of light sources that are induced by the rank-1 latticestratification (reference equation (10)). For each of the “r” rays thatare shot into the scene from the point in the pixel that is associatedwith the “i-th” element ξ_(i) of the random sequence{ξ_(i)}_(i = 0)^(r − 1),if the shot ray intersects an object in the scene, one shadow ray istraced from the point at which the shot ray intersects an object in thescene, to one point within each of the j=1, . . . ,n strata on the lightsource, the point in the stratum that is determined by the “i-th”element of the random sequence {ζ_(i)}_(i = 0)^(r − 1).Since for each element ξ_(i) of the random sequence{ξ_(i)}_(i = 0)^(r − 1)(which, as noted above, is used in determining the location within thepixel from which the ray shot will be shot into the scene), the sameelement ζ_(i) of the random sequence {ζ_(i)}_(i = 0)^(r − 1)will be used in generating values for R_(j)(ζ_(i)) (which is used indetermining the locations in the various strata on the light sourcetowards which the “j” respective shadow rays will be traced) for variousvalues of “j,” the shadow rays that are traced from the point ofincidence of the shot ray will have the same displacements in thevarious strata that are induced on the light source.

On the other hand, for the “r” different points in the pixel from whichrays are shot, locations 10 of which are determined using the differentelements . . . ,ξ_(i−1), ξ_(i), ξ_(i+1) . . . of the random sequence{ξ_(i)}_(i = 0)^(r − 1),a different element . . . ,ζ_(i−1), ζ_(i), ζ_(i+1) . . . of the randomsequence {ζ_(i)}_(i = 0)^(r − 1)will be used in determining the displacements within the various strataon the light source toward which the shadow rays will be traced.Accordingly, the displacements in the respective strata on the lightsource towards which shadow rays will be shot will usually differ forthe different rays that are shot in the scene from the respective pixel.It will be appreciated, however, that, since the elements ζ_(i) are froma random “s₂”-dimensional sequence, it is possible that two moreelements of the sequence {ζ_(i)}_(i = 0)^(r − 1)may have the same or close to the same value, in which case, shadow raysthat are traced from the intersection points in the scene will be tracedto the same or almost the same points on the light source.

The jittered sampling, correlated sampling and trajectory splittingmethodologies as described above will be further described in connectionwith three particular applications that are used in connection withgenerating images in computer graphics, namely, generating photon maps,path tracing and distribution ray tracing.

Photon Map Generation

As noted above, the photon map methodology generally includes twophases, namely, a first “photon map generation” phase followed by a“gather” phase. In the photon map generation phase, the computergraphics system 10 shoots rays that represent paths taken by simulatedphotons emitted by one or more simulated light sources into a simulatedscene. Each rays is in the form of a random walk that is traced throughat least a portion of the simulated scene. The computer graphics systemdetermines whether a ray intersects a surface of an object in the scene,and, if so, stores information relating to the photon associated withthe ray in a database. The information which may include, for example,the photon's color, its direction of incidence, the location of theintersection point, and other information that will be apparent to thoseskilled in the art, in a database. Depending on a number of criteria,including characteristics of the surface of the simulated object at thepoint at which the photon intersects the object, the number of times (ifany) the photon has been previously reflected, and other criteria thatwill be apparent to those skilled in the art, the computer graphicssystem may terminate the photon at the point at which it intersects theobject, or it may allow the photon to be reflected along another ray. Ifthe computer graphics system determines that the photon is to bereflected, the computer graphics system determines the direction ofreflection using a number of criteria, including, among other things,the angle of incidence and characteristics of the simulated object'ssurface, such as whether the surface is specular, glossy or diffuse. Thecomputer graphics system can make use of random factors in several ways,including, for example, determining whether to terminate a ray at asurface that the ray intersects, or to allow the ray it to be reflected.The computer graphics system may also use random factors to jitter theangle with which a ray is to reflected, and as well as in other ways aswill be apparent to those skilled in the art.

After the computer graphics system has shot a number of rays into thescene, it initiates the second “gather” phase of the photon mapmethodology. In that phase, as described above, the computer graphicssystem 10 shoots simulated rays from a simulated image plane of asimulated camera into the scene. The rays that are shot from the imageplane into the scene during the gather phase essentially correspond tothe directions into the scene that are visible from the respectivepoints on the pixel from which they are shot. If a gather ray intersectsthe surface of an object in the scene, the computer graphics systemmakes use of the information regarding the photons, if any, thatintersected the surface proximate the point at which the ray shot fromthe camera intersected the surface, including the directions ofincidence of the respective photons, the direction of incidence of thegather ray, and the various characteristics of the surface on which thegather ray is incident. The computer graphics system uses thatinformation to determine the brightness and color of the scene at thatpoint in the image.

As noted above, for rays that are shot into the simulated scene, randomfactors can be used in several ways, including determining whether toterminate a ray at a surface that the ray intersects or whether to allowthe ray to be reflected, determining the direction of reflection of thephotons that are shot into the scene during the photon map generationphase, and other ways that will be appreciated by those skilled in theart. The random factors that are used to jitter each ray are collectedinto an “s”-dimensional vector ξ ε[0,1)^(s), where “s” is related to themaximum number of reflections or refractions that a photon may havebefore it is terminated. Problems have arisen in the past, sincestratification of high-dimensional integration domains has not beenavailable. However, stratification using rank-1 lattices provides asolution to this problem. In particular, the computer graphics system 10makes use of a sequence of random s-dimensional vectors{ξ_(j)}_(j = 0)^(r − 1),ξ_(j) ε[0,1)^(s), and uses the “j-th” element ξ_(j) of the sequence ingenerating R_(j)(ξ_(j)) (reference equation (10)) for the “j-th” randomwalk. By jittering the random walk using a random number set, thecomputer graphics system 10 can reduce the likelihood that correlationpatterns will develop in the traces that the computer graphics systemshoots into the scene.

Path Tracing

In path tracing, the computer graphics system 10 evaluates a renderingequation that is of the formL({right arrow over (x)},{right arrow over (ω)})+∫_(A) L({right arrowover (x)}′,{right arrow over (ω)}({right arrow over (x)}′,{right arrowover (x)})ƒ_(r)({right arrow over (x)},{right arrow over (ω)}){rightarrow over (x)}′,{right arrow over (x)}))G({right arrow over (x)},{rightarrow over (x)}ω′)V({right arrow over (x)},{right arrow over(x)}′)dA′  (16)where

(i) L({right arrow over (x)},{right arrow over (ω)}) on the left-handside of equation (16) represents the luminance at point {right arrowover (x)} in direction {right arrow over (ω)};

(ii) L_(e)({right arrow over (x)}, {right arrow over (ω)}) on theright-hand side is an emission term representing the illuminationemitted from point {right arrow over (x)} in the direction {right arrowover (ω)};

(iii) L({right arrow over (x)}′, ({right arrow over (x)}′,{right arrowover (x)})) on the right-hand side represents the luminance at anotherpoint {right arrow over (x)}′ in the scene that is reflected in thedirection ω({right arrow over (x)}′,{right arrow over (x)}) toward point{right arrow over (x)};

(iv) ƒ_(r)({right arrow over (x)},ω({right arrow over (x)}′,{right arrowover (x)})) on the right-hand side is a bidirectional scatteringdistribution function that describes how much of the light incident onfrom the point {right arrow over (x)}′ is reflected, refracted orotherwise scattered in the direction ω({right arrow over (x)}′,{rightarrow over (x)}) toward the point {right arrow over (x)}, and isgenerally the sum of a diffuse component, a glossy component and aspecular component;

(v) G({right arrow over (x)},{right arrow over (x)}′) on the right-handside is a geometric term that is a function of the orientation of thesurfaces that contain points {right arrow over (x)}′ and {right arrowover (x)} relative to one another, specifically, $\begin{matrix}{{G\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = \frac{\cos\quad\theta\quad\cos\quad\theta^{\prime}}{{{\overset{->}{x} - {\overset{->}{x}}^{\prime}}}^{2}}} & (17)\end{matrix}$where θ and θ′ are angles of the line between points {right arrow over(x)} and {right arrow over (x)}′, respectively relative to the normalsof the surfaces that contain points {right arrow over (x)} and {rightarrow over (x)}′, and (vi) V({right arrow over (x)},{right arrow over(x)}′) on the right-hand side of equation (16) is a visibility functionthat equals the value one if the point {right arrow over (x)}′ isvisible from the point {right arrow over (x)} and zero if the point{right arrow over (x)}′ is not visible from the point {right arrow over(x)}.It will be appreciated that point {right arrow over (x)} represents thepoint at which the luminance is to be determined and point {right arrowover (x)}′ represents another point in the scene from which light may beemitted (in the case of a light source) or reflected toward the point{right arrow over (x)}. Essentially equation (16) reflects the fact thatthe luminance at point {right arrow over (x)} in the direction {rightarrow over (ω)} is generally the sum of two components. One component,which is represented by item (ii) above (L_(e)({right arrow over(x)},{right arrow over (ω)})), represents the light (if any) that isemitted from the point {right arrow over (x)} in the direction {rightarrow over (ω)}. This component will be non-zero if the point {rightarrow over (x)} is on a light source.

The second component, which is represented by the integral on the righthand side of equation (16), represents the light (if any) that isdirected toward point {right arrow over (x)} from elsewhere in thescene, including, for example, light sources and other surfaces, andthat is directed to or reflected or otherwise scattered towards thepoint {right arrow over (x)} from other points {right arrow over (x)}′in the scene that are visible from point {right arrow over (x)}. Theintegral may be non-zero for points on light sources as well as forpoints on other objects in the scene. The integral in equation (16) isevaluated over the area A′ of the surface of a hemisphere that is formedfrom a sphere that is centered on the point {right arrow over (x)} andgenerally bisected by a plane that is tangent to the surface at point{right arrow over (x)}.

The computer graphics system 10 numerically evaluates the integral inequation (16) by shooting a plurality of rays from respective points onthe pixel into the scene. For each ray shot from the pixel, if the pathof the ray intersects the surface of an simulated object in the scene,in addition to shooting shadow rays to light sources as described aboveto determine the extent to the point on the object is illuminated by alight source, the computer graphics system 10 can also determine theextent to which the point on the object is illuminated by lightreflected from elsewhere in the scene. In the latter operation, thecomputer graphics system 10 can extend the ray by tracing a path fromthe point on the object that the path of the ray as shot from the pixelintersected, in a direction that is generally a reflection of the pathof incidence. The computer graphics system can repeat these operationsthrough a plurality of iterations, or until the path of the ray does notintersect an object in the scene. In addition, the computer graphicssystem 10 will perform similar operations for all of the rays that areshot from the pixel, and the combine the results for the respective raysto determine the pixel value for the pixel.

It will be appreciated that the computer graphics system 10 makes use ofa large number of paths in determining the pixel value, it can usecorrelated sampling to determine the locations of the points of thepixel. Since in correlated sampling the same value ξ is used for all ofthe sample points for the respective pixel, the computer graphics system10 can generate the sample points R_(j)(ξ) for the pixel by use of arigid transformation, which is more efficient than, for example,jittered sampling.

Distribution Ray Tracing

The distribution ray tracing methodology extends the path tracingmethodology as described above, by adding trajectory splitting. In thedistribution ray tracing methodology, if a ray shot from the simulatedimage plane into the scene intersects the simulated surface of asimulated object in a scene, instead of using a single ray reflectedfrom the surface, the computer graphics system reflects a plurality ofrays along a like plurality of paths from the intersection point on thesurface. The computer graphics system divides the hemisphere into aplurality of strata, one for each path that is to be reflected from thesurface. The computer graphics system directs a ray through eachstratum, and can make use of trajectory splitting as described above todetermine the specific direction through the respective strata.Accordingly, the computer graphics system will trace a plurality ofpaths from the point at which the ray's path intersects the point on theobject. The computer graphics system can perform trajectory splitting atany point at which a ray's path intersects a point on an object.

The invention provides a number of advantages. In particular, theinvention provides a new and improved system and computer-implementedmethod for evaluating integrals using a Monte Carlo methodology in whichthe integration domain is partitioned into a plurality of strata using astratification methodology in which the number of strata is independentof the number of dimensions comprising the integration domain.Specifically, the invention provides a methodology that makes use ofstratification based on rank-1 lattices. Systems andcomputer-implemented methods according to the invention find utility ina number of applications, including but not limited to computergraphics.

The invention is also beneficial in connection with one problem that canarise in computer graphics, namely, aliasing, which gives rise toartifacts in the form of, for example, jagged edges. A computer graphicssystem 10 in accordance with the invention does not eliminate aliasingartifacts, but it does, however, mask them by noise that can arise fromrandomization. To accomplish this, the sample points are preferablyrandom as between different pixels; however, the samples can becorrelated within a pixel. The estimator described above in connectionwith the last line of equation (13) matches these requirements. Within apixel, the samples are a shifted lattice, with the amount of shift beingdetermined by the value of the random element ξ that is used for thepixel, providing good uniformity of distribution throughout the area ofthe pixel, which, in turn, provides for fast convergence. For differentpixels, the various values of elements . . . ξ_(i−1),ξ_(i),ξ_(i+1) . . .that are used provides noise that can mask the aliasing artifacts in theimage that is generated.

It will be appreciated that a number of changes and modifications may bemade to the invention as described herein. For example, although theinvention has been described in the context of a computer graphicssystem, it will be appreciated that the invention will find utility inother applications in which integrals are numerically evaluated usingthe Monte Carlo methodology.

In addition, although the invention has been described in connectionwith a rank-1 lattice in the form of a lattice based on the Fibonaccisequence, it will be appreciated that other forms of rank-1 lattices maybe used.

It will be appreciated that a system in accordance with the inventioncan be constructed in whole or in part from special purpose hardware ora general purpose computer system, or any combination thereof, anyportion of which may be controlled by a suitable program. Any programmay in whole or in part comprise part of or be stored on the system in aconventional manner, or it may in whole or in part be provided in to thesystem over a network or other mechanism for transferring information ina conventional manner. In addition, it will be appreciated that thesystem may be operated and/or otherwise controlled by means ofinformation provided by an operator using operator input elements (notshown) which may be connected directly to the system or which maytransfer the information to the system over a network or other mechanismfor transferring information in a conventional manner.

The foregoing description has been limited to a specific embodiment ofthis invention. It will be apparent, however, that various variationsand modifications may be made to the invention, with the attainment ofsome or all of the advantages of the invention. It is the object of theappended claims to cover these and such other variations andmodifications as come within the true spirit and scope of the invention.

1. A system for numerically evaluating an integral of a function over ans-dimensional integration domain comprising: A. a sample point generatorconfigured to generate a selected number of sample points over theintegration domain, the sample points being generated such that there isat least one sample point in each of a plurality of strata distributedover the integration domain, the strata being defined by a rank-1lattice; B. a function value generator configured to, for respectiveones of the sample points, generate a value for the function at therespective sample point; and C. an integral value estimate generatorconfigured to use the function values generated by the function valuegenerator at the respective sample points in generating an estimate forthe value of the integral:
 2. A system as defined in claim 1 in whichthe integration domain comprises the “s”-dimensional unit cube[0,1)^(s).
 3. A system as defined in claim 2 in which, for each stratumA_(j), the sample point generator is configured to generate the “j-th”sample point R_(j)(x) in accordance with the bijectionR_(j) : [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {B\quad x}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 4. A system asdefined in claim 3 in which the sample point generator is configured togenerate the “j-th” sample point such that “x” is a vector.
 5. A systemas defined in claim 4 in which the sample point generator is configuredto generate the “j-th” sample point such that the vector “x” is the samefor all sample points.
 6. A system as defined in claim 4 in which, foreach value of index “j ,” the sample point generator is configured togenerate the “j-th” sample point such that the vector “x” is the “j-th”element of a sequence of vectors P_(n)=(ξ₀, . . . ,ξ_(n-1)).
 7. A systemas defined in claim 6 in which the elements of the sequence aregenerated using a random sequence generation methodology.
 8. A system asdefined in claim 6 in which at least one vector in the sequence is an“s”-dimensional vector.
 9. A system as defined in claim 2 in which thesample point generator is configured to generate the “j-th” sample pointsuch that the “s-”dimensional unit cube [0,1)^(s) comprises an“s₁”-dimensional portion and an “s₂”-dimensional portion, where eachsample point is defined by (ξ_(i),R_(j)(ζ_(i))), where ξ_(i) is anelement of an “s₁”-dimensional sequence, ζ_(i) is an element of an“s₂”-dimensional sequence, and R_(j) is defined by the bijectionR_(j) : [0, 1)^(s₂)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {B\quad x}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 10. A system asdefined in claim 9 in which the sample point generator is configured togenerate the “j-th” sample point such that the elements of the“s₁”-dimensional sequence are generated using a random sequencegeneration methodology.
 11. A system as defined in claim 9 in which thesample point generator is configured to generate the “j-th” sample pointsuch that the elements of the “s₂”-dimensional sequence are generatedusing a random sequence generation methodology.
 12. A computerimplemented method of numerically evaluating an integral of a functionover an s-dimensional integration domain comprising: A. a sample pointgeneration step of generating a selected number of sample points overthe integration domain, the sample points being generated such thatthere is at least one sample point in each of a plurality of stratadistributed over the integration domain, the strata being defined by arank-1 lattice; B. a function value generation step of, for respectiveones of the sample points, generating a value for the function at therespective sample point; and C. an integral value estimate generationstep of using the function values generated during the function valuegeneration step at the respective sample points in generating anestimate for the value of the integral.
 13. A method as defined in claim12 in which the integration domain comprises the “s”-dimensional unitcube [0,1)^(s).
 14. A method as defined in claim 13 in which, for eachstratum A_(j), the sample point generation step includes the step ofgenerating the “j-th” sample point R_(j)(x) in accordance with thebijection R_(j):  [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 15. A method asdefined in claim 14 in which the sample point generation step includesthe step of generating the “j-th” sample point such that “x” is avector.
 16. A method as defined in claim 15 in which the sample pointgeneration step includes the step of generating the “j-th” sample pointsuch that the vector “x” is the same for all sample points.
 17. A methodas defined in claim 15 in which, for each value of index “j,” the samplepoint generation step includes the step of generating the “j-th” samplepoint such that the vector “x” is the “j-th” element of a sequence ofvectors P_(n)=(ξ₀, . . . ,ξ_(n-1)).
 18. A method as defined in claim 17in which the sample point generation step includes the step ofgenerating the elements of the sequence using a random sequencegeneration methodology.
 19. A method as defined in claim 17 in which atleast one vector in the sequence is an “s”-dimensional vector.
 20. Amethod as defined in claim 13 in which the sample point generation stepincludes the step of generating the “j-th” sample point such that the“s-”dimensional unit cube [0,1)^(s) comprises an “s₁”-dimensionalportion and an “s₂”-dimensional portion, where each sample point isdefined by (ξ_(i),R_(j)(ζ_(i))),where ξ_(i) is an element of an“s₁”-dimensional sequence, ζ_(i) is an element of an “s₂”-dimensionalsequence, and R_(j) is defined by the bijection $\begin{matrix}{{R_{j}{\text{:}\quad\left\lbrack {0,1} \right)}^{s_{2}}}->A_{j}} \\\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.\end{matrix}$ where “n” is the selected number of sample points, g is agenerator vector for the rank-1 lattice, and B is a matrix that definesa linear transformation that represents a rotation and/or shear.
 21. Amethod as defined in claim 20 in which the sample point generation stepincludes the step of generating the “j-th” sample point such that theelements of the “s₁”-dimensional sequence are generated using a randomsequence generation methodology.
 22. A method as defined in claim 20 inwhich the sample point generation step includes the step of generatingthe “j-th” sample point such that the elements of the “s₂”-dimensionalsequence are generated using a random sequence generation methodology.23. A computer program product for use with a computer to provide asystem for numerically evaluating an integral of a function over ans-dimensional integration domain, the computer program productcomprising a computer readable medium having encoded thereon: A. asample point generator module configured to enable the computer togenerate a selected number of sample points over the integration domain,the sample points being generated such that there is at least one samplepoint in each of a plurality of strata distributed over the integrationdomain, the strata being defined by a rank-1 lattice; B. a functionvalue generator module configured to enable the computer to, forrespective ones of the sample points, generate a value for the functionat the respective sample point; and C. an integral value estimategenerator module configured to enable the computer to use the functionvalues generated by the function value generator at the respectivesample points in generating an estimate for the value of the integral.24. A computer program product as defined in claim 23 in which theintegration domain comprises the “s”-dimensional unit cube [0,1)^(s).25. A computer program product as defined in claim 24 in which, for eachstratum A_(j), the sample point generator module is configured to enablethe computer to generate the “j-th” sample point R_(j)(x) in accordancewith the bijection R_(j):  [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 26. A computerprogram product as defined in claim 25 in which the sample pointgenerator module is configured to enable the computer to generate the“j-th” sample point such that “x” is a vector.
 27. A computer programproduct as defined in claim 26 in which the sample point generatormodule is configured to enable the computer to generate the “j-th”sample point such that the vector “x” is the same for all sample points.28. A computer program product as defined in claim 26 in which, for eachvalue of index “j,” the sample point generator module is configured toenable the computer to generate the “j-th” sample point such that thevector “x” is the “j-th” element of a sequence of vectors P_(n)=(ξ₀, . .. ,ξ_(n-)1 ).
 29. A computer program product as defined in claim 28 inwhich the elements of the sequence are generated using a random sequencegeneration methodology.
 30. A computer program product as defined inclaim 28 in which at least one vector in the sequence is an“s”-dimensional vector.
 31. A computer program product as defined inclaim 24 in which the sample point generator module is configured toenable the computer to generate the “j-th” sample point such that the“s-”dimensional unit cube [0,1)^(s) comprises an “s₁”-dimensionalportion and an “s₂”-dimensional portion, where each sample point isdefined by (ξ_(i),R_(j)(ζ_(i))), where ξ_(i) is an element of an“s”-dimensional sequence, ζ_(i) is an element of an “s₂”-dimensionalsequence, and R_(j) is defined by the bijection $\begin{matrix}{{R_{j}{\text{:}\quad\left\lbrack {0,1} \right)}^{s_{2}}}->A_{j}} \\\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.\end{matrix}$ where “n” is the selected number of sample points, g is agenerator vector for the rank-1 lattice, and B is a matrix that definesa linear transformation that represents a rotation and/or shear.
 32. Acomputer program product as defined in claim 31 in which the samplepoint generator module is configured to enable the computer to generatethe “j-th” sample point such that the elements of the “s₁”-dimensionalsequence are generated using a random sequence generation methodology.33. A computer program product as defined in claim 31 in which thesample point generator module is configured to enable the computer togenerate the “j-th” sample point such that the elements of the“s₂”-dimensional sequence are generated using a random sequencegeneration methodology.
 34. A computer graphics system for generating apixel value for a pixel in an image, the pixel being representative of apoint in a scene, the scene comprising at least one object and at leastone light source, the computer graphics system generating the pixelvalue by an evaluation of an integral of a selected function over ans-dimensional integration domain comprising: A. a sample point generatorconfigured to generate a selected number of sample points over theintegration domain, the sample points being generated such that there isat least one sample point in each of a plurality of strata distributedover the integration domain, the strata being defined by a rank-1lattice; B. a function value generator configured to, for respectiveones of the sample points, generate a value for the function at therespective sample point; and C. an integral value estimate generatorconfigured to use the function values generated by the function valuegenerator at the respective sample points in generating an estimate forthe value of the integral in relation to the at least one object and atleast one light source, the estimate corresponding to the pixel valuefor the image.
 35. A system as defined in claim 34 in which theintegration domain comprises the “s”-dimensional unit cube [0,1)^(s) andfor each stratum A_(j), the sample point generator is configured togenerate the “j-th” sample point R_(j)(x) in accordance with thebijection R_(j):  [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 36. A system asdefined in claim 35 in which the sample point generator is configured touse the same vector “x” for at least some of the sample points.
 37. Asystem as defined in claim 36 in which at least part of the selectedfunction represents a path tracing methodology in which the functionvalue generator is configured to generate the function value associatedwith a point in the scene in relation to a plurality of simulated tracesprojected from the point, the directions of the traces being determinedby strata defined on at least a portion of a sphere centered on thepoint induced by at least a portion of the rank-1 lattice.
 38. A systemas defined in claim 35 in which, for each value of index “j,” the samplepoint generator is configured to generate the “j-th” sample point suchthat the vector “x” is the “j-th” element of a sequence of vectorsP_(n)=(ξ₀, . . . ,ξ_(n-1)).
 39. A system as defined in claim 38 in whichthe elements of the sequence are generated using a random sequencegeneration methodology.
 40. A system as defined in claim 39 in which atleast part of the selected function represents a photon map methodologyin which the function value generator is configured to generate thefunction values by using sample points generated by the sample pointgenerator module using least some of the vectors from the sequence tocontrol a random walk of traces from the light source through at least aportion of the scene.
 41. A system as defined in claim 35 in which thesample point generator is configured to generate the “j-th” sample pointsuch that the “s-” dimensional unit cube [0,1)^(s) comprises an“s₁”-dimensional portion and an “s₂”-dimensional portion, where eachsample point is defined by (ξ_(i),R_(j)(ζ_(i))), where ξ_(i) is anelement of an “s₁”-dimensional sequence, ζ_(i) is an element of an“s₂”-dimensional sequence, and R_(j) is defined by the bijection$\begin{matrix}{{R_{j}{\text{:}\quad\left\lbrack {0,1} \right)}^{s_{2}}}->A_{j}} \\\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.\end{matrix}$ where “n” is the selected number of sample points, g is agenerator vector for the rank-1 lattice, and B is a matrix that definesa linear transformation that represents a rotation and/or shear.
 42. Asystem as defined in claim 41 in which the sample point generator isconfigured to generate the “j-th” sample point such that the elements ofthe “s₁”-dimensional sequence are generated using a random sequencegeneration methodology.
 43. A system as defined in claim 41 in which thesample point generator is configured to generate the “j-th” sample pointsuch that the elements of the “s₂”-dimensional sequence are generatedusing a random sequence generation methodology.
 44. A system as definedin claim 41 in which at least part of the selected function represents ashadow ray methodology, in which a plurality of shadow rays comprisingsimulated traces are traced from a point in the scene towards respectivepoints on the light source, the locations of the points on the lightsource being determined by R_(j), the function value generator beingconfigured to generate the function value in relation to whetherrespective ones of the shadow rays intersect at least one of saidobjects in the scene.
 45. A system as defined in claim 41 in which atleast a part of the selected function represents a distribution pathtracing methodology, in which the function value generator is configuredto generate the function value associated with a point in the scene inrelation to a plurality of simulated traces projected from the point,the directions of the traces being determined by strata defined on atleast a portion of a sphere centered on the point induced by R_(j). 47.A computer implemented method of generating a pixel value for a pixel inan image, the pixel being representative of a point in a scene, thescene comprising at least one object and at least one light source, thecomputer graphics method generating the pixel value by an evaluation ofan integral of a selected function over an s-dimensional integrationdomain, the computer graphics method comprising: A. a sample pointgeneration step of generating a selected number of sample points overthe integration domain, the sample points being generated such thatthere is at least one sample point in each of a plurality of stratadistributed over the integration domain, the strata being defined by arank-1 lattice; B. a function value generation step of, for respectiveones of the sample points, generating a value for the function at therespective sample point; and C. an integral value estimate generationstep of using the function values generated during the function valuegeneration step at the respective sample points in generating anestimate for the value of the integral in relation to the at least oneobject and at least one light source, the estimate corresponding to thepixel value for the image.
 48. A method as defined in claim 47 in whichthe integration domain comprises the “s”-dimensional unit cube [0,1)^(s)and for each stratum A_(j), the sample point generation step includingthe step of generating the “j-th” sample point R_(j)(x) in accordancewith the bijection R_(j):  [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 49. A method asdefined in claim 48 in which the sample point generation step includesthe step of using the same vector “x” for at least some of the samplepoints.
 50. A method as defined in claim 49 in which at least part ofthe selected function represents a path tracing methodology in which thefunction value generation step includes the step of generating thefunction value associated with a point in the scene in relation to aplurality of simulated traces projected from the point, the directionsof the traces being determined by strata defined on at least a portionof a sphere centered on the point induced by at least a portion of therank-1 lattice.
 51. A method as defined in claim 48 in which, for eachvalue of index “j,” the sample point generation step includes the stepof generating the “j-th” sample point such that the vector “x” is the“j-th” element of a sequence of vectors P_(n)=(ξ₀, . . .ξ_(n-1)).
 52. Amethod as defined in claim 51 in which the sample point generation stepincludes the step of generating the elements of the sequence using arandom sequence generation methodology.
 53. A method as defined in claim52 in which at least part of the selected function represents a photonmap methodology in which the function value generation step includes thestep of generating the function values by using sample points generatedduring the sample point generation step using least some of the vectorsfrom the sequence to control a random walk of traces from the lightsource through at least a portion of the scene.
 54. A method as definedin claim 48 in which the sample point generation step includes the stepof generating the “j-th” sample point such that the “s-”dimensional unitcube [0,1)^(s) comprises an “s₁”-dimensional portion and an“s₂”-dimensional portion, where each sample point is defined by(ξ_(i),R_(j)(ζ_(i))), where ζ_(i) is an element of an “s₁”-dimensionalsequence, ζ_(i) is an element of an ”s₂”-dimensional sequence, and R_(j)is defined by the bijection $\begin{matrix}{{R_{j}{\text{:}\quad\left\lbrack {0,1} \right)}^{s_{2}}}->A_{j}} \\\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.\end{matrix}$ where “n” is the selected number of sample points, g is agenerator vector for the rank-1 lattice, and B is a matrix that definesa linear transformation that represents a rotation and/or shear.
 55. Amethod as defined in claim 54 in which the sample point generation stepincludes the step of generating the “j-th” sample point such that theelements of the “s₁”-dimensional sequence are generated using a randomsequence generation methodology.
 56. A method as defined in claim 54 inwhich the sample point generation step includes the step of generatingthe “j-th” sample point such that the elements of the “s₂”-dimensionalsequence are generated using a random sequence generation methodology.57. A method as defined in claim 54 in which at least part of theselected function represents a shadow ray methodology, in which aplurality of shadow rays comprising simulated traces are traced from apoint in the scene towards respective points on the light source, thelocations of the points on the light source being determined by R_(j),the function value generation step including the step of generating thefunction value in relation to whether respective ones of the shadow raysintersect at least one of said objects in the scene.
 58. A method asdefined in claim 54 in which at least a part of the selected functionrepresents a distribution path tracing methodology, in which thefunction value generation step includes the step of generating thefunction value associated with a point in the scene in relation to aplurality of simulated traces projected from the point, the directionsof the traces being determined by strata defined on at least a portionof a sphere centered on the point induced by R_(j).
 59. A computerprogram product for use in connection with a computer to provide acomputer graphics system for generating a pixel value for a pixel in animage, the pixel being representative of a point in a scene, the scenecomprising at least one object and at least one light source, thecomputer graphics system generating the pixel value by an evaluation ofan integral of a selected function over an s-dimensional integrationdomain, the computer program product comprising a computer-readablemedium having encoded hereon: A. a sample point generator moduleconfigured to enable the computer to generate a selected number ofsample points over the integration domain, the sample points beinggenerated such that there is at least one sample point in each of aplurality of strata distributed over the integration domain, the stratabeing defined by a rank-1 lattice; B. a function value generator moduleconfigured to enable the computer to, for respective ones of the samplepoints, generate a value for the function at the respective samplepoint; and C. an integral value estimate generator module configured toenable the computer to use the function values generated by the functionvalue generator at the respective sample points in generating anestimate for the value of the integral in relation to the at least oneobject and at least one light source, the estimate corresponding to thepixel value for the image.
 60. A computer program product as defined inclaim 59 in which the integration domain comprises the “s”-dimensionalunit cube [0,1)^(s) and for each stratum A_(j), the sample pointgenerator module is configured to enable the computer to generate the“j-th” sample point R_(j)(x) in accordance with the bijectionR_(j):  [0, 1)^(s)− > A_(j)$\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s}} \right.$where “n” is the selected number of sample points, g is a generatorvector for the rank-1 lattice, and B is a matrix that defines a lineartransformation that represents a rotation and/or shear.
 61. A computerprogram product as defined in claim 60 in which the sample pointgenerator module is configured to enable the computer to use the samevector “x” for at least some of the sample points.
 62. A computerprogram product as defined in claim 61 in which at least part of theselected function represents a path tracing methodology in which thefunction value generator module is configured to enable the computer togenerate the function value associated with a point in the scene inrelation to a plurality of simulated traces projected from the point,the directions of the traces being determined by strata defined on atleast a portion of a sphere centered on the point induced by at least aportion of the rank-1 lattice.
 63. A computer program product as definedin claim 60 in which, for each value of index “j,” the sample pointgenerator module is configured to enable the computer to generate the“j-th” sample point such that the vector “x” is the “j-th” element of asequence of vectors P_(n)=(ξ₎, . . . ,ξ_(n-1)).
 64. A computer programproduct as defined in claim 63 in which the elements of the sequence aregenerated using a random sequence generation methodology.
 65. A computerprogram product as defined in claim 64 in which at least part of theselected function represents a photon map methodology in which thefunction value generator module is configured to enable the computer togenerate the function values by using sample points generated by thesample point generator module using least some of the vectors from thesequence to control a random walk of traces from the light sourcethrough at least a portion of the scene.
 66. A computer program productas defined in claim 60 in which the sample point generator module isconfigured to enable the computer to generate the “j-th” sample pointsuch that the “s-”dimensional unit cube [0,1)^(s) comprises an“s₁”-dimensional portion and an “s₂”-dimensional portion, where eachsample point is defined by (ξ_(i),R_(j)(ζ_(i))), where ζ_(i) is anelement of an “s₁”-dimensional sequence, ζ_(i) is an element of an“s₂”-dimensional sequence, and R_(j) is defined by the bijection$\begin{matrix}{{R_{j}{\text{:}\quad\left\lbrack {0,1} \right)}^{s_{2}}}->A_{j}} \\\left. x\mapsto{\left( {{\frac{j}{n} \cdot g} + {Bx}} \right)\quad{{mod}\quad\left\lbrack {0,1} \right)}^{s_{2}}} \right.\end{matrix}$ where “n” is the selected number of sample points, g is agenerator vector for the rank-1 lattice, and B is a matrix that definesa linear transformation that represents a rotation and/or shear.
 67. Acomputer program product as defined in claim 66 in which the samplepoint generator module is configured to enable the computer to generatethe “j-th” sample point such that the elements of the “s₁”-dimensionalsequence are generated using a random sequence generation methodology.68. A computer program product as defined in claim 66 in which thesample point generator module is configured to enable the computer togenerate the “j-th” sample point such that the elements of the“s₂”-dimensional sequence are generated using a random sequencegeneration methodology.
 69. A computer program product as defined inclaim 66 in which at least part of the selected function represents ashadow ray methodology, in which a plurality of shadow rays comprisingsimulated traces are traced from a point in the scene towards respectivepoints on the light source, the locations of the points on the lightsource being determined by R_(j), the function value generator beingconfigured to generate the function value in relation to whetherrespective ones of the shadow rays intersect at least one of saidobjects in the scene.
 70. A computer program product as defined in claim66 in which at least a part of the selected function represents adistribution path tracing methodology, in which the function valuegenerator module is configured to enable the computer to generate thefunction value associated with a point in the scene in relation to aplurality of simulated traces projected from the point, the directionsof the traces being determined by strata defined on at least a portionof a sphere centered on the point induced by R_(j).